file1=read.table("E:/2_5hmc_yjp_bam/5hmc_file/2_5hmc_yjp_bam/DMR_edgeR_noRMS/6v6p_8v8p/chipseq_encode/motif_overlap_TF_chipseq.txt",head=T,sep="\t")
#file1=file1[file1$BF_in_DC>10,]
ID=as.character(unique(file1$unitID))

sel=list.files(pattern="bayes_p.txt")
sel1=c("X2B_X1T","M8_M7","M6_M5","M2_M1","M48_M47","M49_M50")
for (i in 3:length(sel)){
file=read.table(sel[i],head=T,sep="\t")
file=file[file$exits=="TRUE",]
file$unitID=paste0(file$chrom,":",file$position)
file=file[file$unitID %in%ID,]
file$tumor_var_freq=as.numeric(gsub("%","",file$tumor_var_freq))/100
file$normal_var_freq=as.numeric(gsub("%","",file$normal_var_freq))/100
file$total_reads=as.numeric(file$normal_reads1)+as.numeric(file$normal_reads2)+as.numeric(file$tumor_reads1)+as.numeric(file$tumor_reads2)
file=file[file$total_reads>=10,]
#file$DU=abs(file$normal_var_freq-0.5)
#file$DA=abs(file$tumor_var_freq-0.5)
file$p.diff=lapply(1:nrow(file),function(x) fisher.test(rbind(c(file[x,5],file[x,6]),c(file[x,9],file[x,10])))$p.value)#计算各种统计值
#file$p.biasA=lapply(1:nrow(file),  function(x) binom.test(c(file[x,9],file[x,10]),p = 0.5)$p.value)
#file$p.biasU=lapply(1:nrow(file),  function(x) binom.test(c(file[x,5],file[x,6]),p = 0.5)$p.value)
tmp=data.frame(file[,5:7],file[,9:11],p.diff=as.numeric(file$p.diff),unitID=file$unitID)
fn=unlist(strsplit(gsub(".snp.bayes_p.txt","",sel[i]),"_"))
names(tmp)=c(paste0(fn[1],c("_reads1","_reads2","_var_freq")),paste0(fn[2],c("_reads1","_reads2","_var_freq")),paste0(fn[1],"_",fn[2],"_p.diff"),"unitID")
final=merge(final,tmp,by="unitID",all=T)
}



library(metaSeq)
file=read.table("AShM_for_meta.txt",head=T,sep="\t")
allsample=gsub("_reads2","",colnames(file)[grep("reads2",colnames(file))])
sel1=c("X2B_X1T","M8_M7","M6_M5","M2_M1","M48_M47","M50_M49","M18_M17","M20_M19","M22_M21","M40_M39")###单发和健康的
#sel1=c("X2B_X1T","M8_M7","M50_M49","M52_M51")
#sel1=c("M6_M5","M2_M1","M42_M41","M44_M43","M48_M47")
#sel1=c("X2B_X1T","M8_M7","M6_M5","M2_M1","M48_M47","M50_M49","M30_M29","M26_M25","M35_M36","M28_M27")  ####包括双发
sel1=c("M30_M29","M26_M25","M35_M36","M28_M27","M18_M17","M20_M19","M22_M21","M40_M39")
sel1=unlist(strsplit(sel1,"_"))

for (i in 1:nrow(file)){
mta1=as.numeric(file[i,paste(sel1[seq(2,length(sel1),2)],"_var_freq",sep="")])
mta2=as.numeric(file[i,paste(sel1[seq(1,length(sel1),2)],"_var_freq",sep="")])
mtp=as.numeric(file[i,paste(sel1[seq(1,length(sel1),2)],"_",sel1[seq(2,length(sel1),2)],"_p.diff",sep="")])
mta1=mta1[!is.na(mta1)]
mta2=mta2[!is.na(mta2)]
mtp=mtp[!is.na(mtp)]
if (length(mtp)>0){
mtp1=mtp/2
mtp2=1-mtp1
pair_up=mtp1
pair_down=mtp2
if(length(which(mta1<=mta2)) >= 1){
pair_up[which(mta1<=mta2)]=mtp2[which(mta1<=mta2)]
pair_down[which(mta1<=mta2)]=mtp1[which(mta1<=mta2)]
}
upper=matrix(pair_up,nrow=1)
lower=matrix(pair_down,nrow=1)
weight=rep(1,length(mtp))
if (length(weight)>=2){
result2=other.oneside.pvalues(upper,lower,weight)
S=Stouffer.test(result2)
metap=S$Upper
metap[which(S$Upper>S$Lower)]=S$Lower[which(S$Upper>S$Lower)]
file[i,"metap_in_DC"]=metap
}
if (length(weight)==1){
file[i,"metap_in_DC"]=mtp
}
}
}
file$metap_in_DC_fdr=p.adjust(file$metap_in_DC,method = "BH")
write.table(file,"ASM_meta.txt",quote=F,row.names=F,sep="\t")
